Chapter 11 Dominance-microbiome analysis

load("data/data.Rdata")
load("hmsc/fit_model1_250_10.Rdata")

11.0.1 Posterior estimates

The next step is to get the posterior estimates of beta, which is the parameter that links dominance with the microbiome. We employ a support threshold of ‘0.05’ to assign significance, meaning that parameters with posterior overlaps of <10% are considered significantly different.

# Select desired support threshold
support=0.9
negsupport=1-support

# Basal tree
postestimates_tree <- genome_tree %>%
        keep.tip(., tip=m$spNames) 

# Posterior estimate table
post_beta <- getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(genome=factor(genome, levels=rev(postestimates_tree$tip.label))) %>%
    mutate(value = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    mutate(value=factor(value, levels=c("Positive","Neutral","Negative"))) %>%
    pivot_wider(names_from = variable, values_from = value) %>%
    rename(intercept=2) %>%
    select(genome,dominance,groupvariable) %>%
    column_to_rownames(var="genome")

# Aggregate basal GIFT into elements
function_table <- genome_gifts %>%
    to.elements(., GIFT_db) %>%
    as.data.frame()

#Phylums
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% postestimates_tree$tip.label) %>%
    arrange(match(genome, postestimates_tree$tip.label)) %>%
    mutate(phylum = factor(phylum, levels = unique(phylum))) %>%
    column_to_rownames(var = "genome") %>%
    select(phylum)


colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
  right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% postestimates_tree$tip.label) %>%
    arrange(match(genome, postestimates_tree$tip.label)) %>%
     select(phylum, colors) %>%
    unique() %>%
    arrange(phylum) %>%
    select(colors) %>%
    pull()

# Basal ggtree
postestimates_tree <- postestimates_tree %>%
        force.ultrametric(.,method="extend") %>%
        ggtree(., size = 0.3)

***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************

#Add phylum colors next to the tree tips
postestimates_tree <- gheatmap(postestimates_tree, phylum_colors, offset=-0.2, width=0.1, colnames=FALSE) +
      scale_fill_manual(values=colors_alphabetic)+
      labs(fill="Phylum")

#Reset fill scale to use a different colour profile in the heatmap
postestimates_tree <- postestimates_tree + new_scale_fill()

# Add posterior significant heatmap
postestimates_tree <- gheatmap(postestimates_tree, post_beta, offset=0, width=0.2, colnames=TRUE, colnames_position="top",colnames_angle=90, colnames_offset_y=1, hjust=0) +
        scale_fill_manual(values=c("#be3e2b","#f4f4f4","#b2b530"))+
        labs(fill="Trend")

#Reset fill scale to use a different colour profile in the heatmap
postestimates_tree <- postestimates_tree + new_scale_fill()

#Add functions heatmap
postestimates_tree <- gheatmap(postestimates_tree, function_table, offset=0.6, width=3.5, colnames=FALSE) +
            vexpand(.08) +
            coord_cartesian(clip = "off") +
            scale_fill_gradient(low = "#f4f4f4", high = "steelblue", na.value="white") +
            labs(fill="Function fullness")

#Reset fill scale to use a different colour profile in the heatmap
postestimates_tree <- postestimates_tree + new_scale_fill()

# Add completeness barplots
postestimates_tree <- postestimates_tree +
            geom_fruit(data=genome_metadata,
            geom=geom_bar,
            grid.params=list(axis="x", text.size=2, nbreak = 1),
            axis.params=list(vline=TRUE),
            mapping = aes(x=length, y=genome, fill=completeness),
                 offset = 3.85,
                 orientation="y",
                 stat="identity") +
            scale_fill_gradient(low = "#cf8888", high = "#a2cc87") +
            labs(fill="Genome completeness")

postestimates_tree +
        vexpand(.25, 1) # expand top 

11.0.2 Positively associated genomes

getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(trend = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    filter(variable=="dominance", trend=="Positive") %>%
    arrange(-value) %>%
    left_join(genome_metadata,by=join_by(genome==genome)) %>%
    select(genome,phylum,class,order,species,value) %>%
    tt()

genome	phylum	class	order	species	value
A_DRC05_bin.59	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__Duncaniella sp910577345	1.000
C_HT_bin.102	p__Bacillota	c__Bacilli	o__Lactobacillales	s__Limosilactobacillus reuteri	0.998
A_HTC12_bin.48	p__Bacillota	c__Bacilli	o__Lactobacillales	s__Lactobacillus johnsonii	0.995
C_HT_bin.9	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__QXXE01 sp910589115	0.994
A_DTC03_bin.77	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__CAG-873 sp910586975	0.992
B_C04M3_bin.80	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__	0.991
B_C05M2_bin.12	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__CAG-873 sp910578095	0.991
A_CDC10_bin.30	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__Duncaniella sp910589485	0.990
B_C03M5_bin.33	p__Bacillota_B	c__Dehalobacteriia	o__UBA4068	s__	0.990
A_CDC13_bin.65	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Eubacterium_J sp009774535	0.987
A_T3C12_bin.34	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__CAG-485 sp009775355	0.985
C_DR_bin.56	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__CAG-485 sp910578075	0.985
B_C04M3_bin.108	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Acutalibacter sp009936035	0.984
B_C06M3_bin.61	p__Bacillota	c__Bacilli	o__Erysipelotrichales	s__C-19 sp910576855	0.981
B_C12F1_bin.72	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__Duncaniella sp910578515	0.980
C_HR_bin.153	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__	0.980
A_T1C10_bin.3	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__UBA7173 sp002491305	0.979
B_C11F5_bin.68	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__UBA7173 sp001689485	0.979
B_C12F1_bin.61	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__CAG-873 sp009775535	0.977
A_DRC04_bin.49	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Enterenecus sp910587255	0.975
C_DR_bin.4	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Avispirillum sp011957885	0.974
A_HTC03_bin.32	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Pelethomonas sp910587645	0.973
B_C04M5_bin.101	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__Paramuribaculum sp910579675	0.973
B_C11F3_bin.29	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__CAG-95 sp009911035	0.973
A_CDC12_bin.65	p__Bacillota_A	c__Clostridia	o__Christensenellales	s__Coproplasma sp003979335	0.972
B_C13F3_bin.77	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__CAG-873 sp910587915	0.972
A_CDC13_bin.13	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__UBA3402 sp002358555	0.971
B_C05M3_bin.38	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Eubacterium_J sp910574915	0.971
B_C13F3_bin.53	p__Bacillota	c__Bacilli	o__Erysipelotrichales	s__Erysipelatoclostridium cocleatum	0.970
A_DTC13_bin.55	p__Bacillota_A	c__Clostridia	o__Christensenellales	s__RACS-045 sp910579785	0.969
A_T2C12_bin.11	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__NSJ-51 sp910585415	0.969
B_C06M5_bin.62	p__Bacillota_A	c__Clostridia	o__Christensenellales	s__CAG-475 sp910577815	0.968
B_C13F3_bin.5	p__Bacillota_A	c__Clostridia	o__TANB77	s__	0.967
B_C13F3_bin.97	p__Bacillota_A	c__Clostridia	o__TANB77	s__	0.967
C_CR_bin.18	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__VSOB01 sp910589075	0.966
C_DT_bin.159	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__	0.965
A_T2C05_bin.3	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__Muribaculum gordoncarteri	0.964
B_C13F2_bin.43	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Enterenecus sp910587255	0.963
B_C04M3_bin.128	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Anaerotruncus sp000403395	0.962
B_C13F3_bin.35	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__MGBC136627 sp910585935	0.961
A_OPC13_bin.50	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__UBA1405 sp910589145	0.960
B_C05M1_bin.75	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Clostridium_Q sp009911305	0.959
B_C13F3_bin.70	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__	0.959
A_CDC05_bin.89	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__	0.958
A_OPC10_bin.79	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__UBA3402 sp910574345	0.957
C_OP_bin.44	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Ventrimonas sp910586205	0.957
C_HT_bin.137	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__	0.956
C_CR_bin.157	p__Bacillota_A	c__Clostridia	o__Christensenellales	s__	0.954
C_T3_bin.64	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__NSJ-51 sp910585415	0.954
C_DR_bin.31	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Anaerotignum sp910576545	0.953
A_OPC10_bin.40	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__	0.952
B_C12F2_bin.47	p__Bacillota	c__Bacilli	o__Lactobacillales	s__Enterococcus_D casseliflavus	0.952
C_HR_bin.13	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Evtepia sp910584805	0.952
C_T2_bin.55	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__	0.952
B_C04M2_bin.33	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Enterenecus sp910585265	0.951
B_C04M3_bin.31	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Angelakisella sp013316495	0.951
A_OPC03_bin.64	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Acetatifactor sp910589655	0.950
A_T1C11_bin.9	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__CAG-485 sp002362485	0.949
B_C04M3_bin.18	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Eubacterium_J sp910586445	0.949
B_C03M3_bin.127	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Acutalibacter sp009936035	0.948
B_C04M4_bin.45	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__	0.948
A_T3C03_bin.17	p__Bacillota_A	c__Clostridia	o__Christensenellales	s__QANA01 sp910588425	0.947
B_C04M2_bin.76	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__CAG-56 sp004793585	0.945
A_HRC11_bin.35	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__	0.944
B_C04M4_bin.5	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Choladocola sp009774145	0.944
C_DT_bin.37	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Acutalibacter sp009936035	0.944
C_DR_bin.21	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__MD308 sp910578455	0.943
A_CDC03_bin.78	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Lawsonibacter sp910577805	0.942
B_C13F5_bin.28	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__Duncaniella sp910576785	0.942
A_HRC10_bin.28	p__Bacillota_A	c__Clostridia	o__Christensenellales	s__Gallimonas sp910578065	0.941
A_T3C11_bin.23	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__RGIG4284 sp910576205	0.941
A_HTC12_bin.21	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__UBA3402 sp910575585	0.940
A_T2C05_bin.28	p__Bacillota	c__Bacilli	o__Erysipelotrichales	s__MGBC163490 sp910588075	0.940
B_C10F2_bin.5	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__	0.940
C_HR_bin.146	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Choladocola sp910575445	0.940
B_C05M5_bin.57	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Pelethomonas sp910587785	0.939
B_C10F5_bin.45	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Dysosmobacter sp000403435	0.939
B_C12F3_bin.66	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__	0.939
B_C11F5_bin.87	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__	0.938
A_CRC05_bin.32	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Eubacterium_J sp910579075	0.937
B_C04M1_bin.34	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Pelethomonas sp910579965	0.937
B_C04M3_bin.90	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__14-2 sp910580015	0.937
C_T2_bin.30	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__	0.934
C_T3_bin.41	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Acetatifactor sp910584235	0.934
C_CD_bin.124	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Pelethomonas sp910579965	0.933
C_OP_bin.73	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__UBA3402 sp910585435	0.933
C_HR_bin.50	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__UBA3402 sp910588325	0.932
A_CDC03_bin.64	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__	0.930
B_C03M3_bin.121	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Lachnoclostridium_B sp910574185	0.930
B_C10F1_bin.59	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Choladocola sp910583895	0.929
B_C10F5_bin.36	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__UBA3282 sp003611805	0.928
B_C11F5_bin.21	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__COE1 sp000403335	0.928
A_CDC10_bin.13	p__Bacillota_A	c__Clostridia	o__Christensenellales	s__MGBC113161 sp910587565	0.927
A_DRC13_bin.41	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Acutalibacter muris	0.927
C_CD_bin.62	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__1XD42-69 sp009911505	0.925
C_OP_bin.15	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Roseburia sp910584215	0.925
A_CDC05_bin.91	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__	0.923
C_T2_bin.92	p__Bacillota_A	c__Clostridia	o__TANB77	s__	0.923
A_CDC05_bin.30	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__	0.922
A_CRC04_bin.71	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Sporofaciens sp910574885	0.922
B_C13F3_bin.99	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Merdisoma sp910574255	0.922
B_C04M1_bin.8	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Merdisoma sp910576325	0.921
B_C03M2_bin.44	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__Duncaniella muris	0.920
B_C06M5_bin.16	p__Actinomycetota	c__Coriobacteriia	o__Coriobacteriales	s__Adlercreutzia muris	0.920
C_T3_bin.99	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Pelethomonas sp910577915	0.920
B_C05M1_bin.6	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Ventrimonas sp910584845	0.919
A_DRC06_bin.37	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__CAG-873 sp009775265	0.918
B_C12F1_bin.38	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Dysosmobacter sp910588005	0.918
B_C10F1_bin.10	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Ventrimonas sp910577765	0.917
B_C13F3_bin.93	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Anaerotruncus sp003612625	0.917
A_HRC05_bin.71	p__Bacillota_A	c__Clostridia	o__Christensenellales	s__MGBC100320 sp910577995	0.915
A_T3C13_bin.3	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__	0.914
B_C03M1_bin.60	p__Bacillota_A	c__Clostridia	o__Christensenellales	s__MGBC100798 sp910584975	0.911
C_OP_bin.100	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__	0.910
C_DT_bin.155	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__	0.909
A_T1C03_bin.24	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__UMGS1370 sp910577645	0.907
B_C12F1_bin.23	p__Bacillota_A	c__Clostridia	o__Oscillospirales	s__Pelethomonas sp910579965	0.907
A_T2C13_bin.14	p__Bacteroidota	c__Bacteroidia	o__Bacteroidales	s__UBA7173 sp001689685	0.906
A_CRC05_bin.73	p__Bacillota_A	c__Clostridia	o__Christensenellales	s__MGBC100320 sp910588305	0.905
B_C03M1_bin.29	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Ventrimonas sp009911065	0.905
B_C04M5_bin.14	p__Actinomycetota	c__Coriobacteriia	o__Coriobacteriales	s__Adlercreutzia caecimuris	0.905
C_HR_bin.28	p__Bacillota_A	c__Clostridia	o__Lachnospirales	s__Sporofaciens sp910585725	0.903
B_C04M5_bin.93	p__Bacillota_A	c__Clostridia	o__Christensenellales	s__Gallimonas sp910588465	0.902

11.0.3 Negatively associated genomes

getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(trend = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    filter(variable=="dominance", trend=="Negative") %>%
    arrange(value) %>%
    left_join(genome_metadata,by=join_by(genome==genome)) %>%
    select(genome,phylum,class,order,species,value) %>%
    tt()

genome	phylum	class	order	species	value
A_T3C05_bin.41	p__Pseudomonadota	c__Gammaproteobacteria	o__Burkholderiales	s__Parasutterella excrementihominis	0.095

11.1 Predict responses to dominance

Based on the fitted model, we can predict the structure of the microbial community for each level of dominance for a certain value of sequencing depth.

# Select modelchain of interest
load("hmsc/fit_model1_250_10.Rdata")

gradient = seq(0,1,0.1)
gradientlength = length(gradient)

#Treatment-specific gradient predictions
pred <- constructGradient(m, 
                      focalVariable = "dominance", 
                      non.focalVariables = list(logseqdepth=list(1)), 
                      ngrid=gradientlength) %>%
             predict(m, Gradient = ., expected = TRUE)

# weights:  3 (2 variable)
initial  value 162.196440 
final  value 133.205769 
converged

predY <- pred %>%
        do.call(rbind,.) %>%
        as.data.frame() %>%
        mutate(dominance=rep(gradient,1000)) %>%
        pivot_longer(-c(dominance), names_to = "genome", values_to = "value")

11.1.1 Plot responses to dominance

The estimated responses of genomes exhibiting significant positive and negative association with dominance.

#Get phylum colors from the EHI standard
phylum_colors <- genome_metadata %>%
    left_join(read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv"), by=join_by(phylum == phylum)) %>%
    arrange(match(genome, genome_tree$tip.label)) %>%
    select(phylum, colors) %>%
    unique() %>%
    arrange(phylum) %>%
    #slice(2:5) %>%
    select(colors) %>%
    pull()

getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(trend = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    filter(variable=="dominance") %>%
    select(genome,trend) %>%
    left_join(predY, by=join_by(genome==genome)) %>%
    #filter(trend != "Neutral") %>%
    #filter(genome %in% predY_asc)  %>% #only display genomes with contrasting dynamics across treatments
    group_by(genome, trend, dominance) %>%
    summarize(value = mean(value, na.rm = TRUE)) %>%
    left_join(genome_metadata, by=join_by(genome == genome)) %>%
    ggplot(aes(x=dominance, y=value, group=genome, color=phylum, linetype=trend)) + 
        geom_line() +
        scale_linetype_manual(values=c("solid","dashed","solid")) +
        scale_color_manual(values=phylum_colors) +
        facet_grid(fct_rev(trend) ~ phylum) +
        labs(y="Genome abundance (log)",x="Dominance") +
        theme(legend.position = "none") +
        theme_minimal() +
         theme(legend.position = "none",
               axis.text.x = element_text(angle = 45, hjust = 0.8,),
               axis.line.x = element_line(size = 0.5, linetype = "solid", colour = "black"),
)

11.2 Plot traits linked to trends

11.2.1 Test functional trait differences

Using Wilcoxon test with Benjamini-Hochberg procedure for FDR, identify functional traits that differ between genomes that are either positively or negatively association with dominance.

postestimate_test <- getPostEstimate(hM=m, parName="Beta")$support %>%
      as.data.frame() %>%
      mutate(variable=m$covNames) %>%
      pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
      mutate(trend = case_when(
            value >= support ~ "Positive",
            value <= negsupport ~ "Negative",
            TRUE ~ "Neutral")) %>%
      filter(variable=="dominance", trend!="Neutral") %>%
      select(genome,trend) %>%
      left_join(function_table %>% rownames_to_column(var="genome"), by=join_by(genome==genome)) %>%
      pivot_longer(-c(genome,trend), names_to = "trait", values_to = "value") %>%
      group_by(trait) %>%
      summarise(p_value = wilcox.test(value ~ trend)$p.value) %>%
      mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
      filter(p_adjust < 0.05)

11.2.2 Community functional predictions

11.2.2.1 Element level

elements_table <- genome_gifts %>%
    to.elements(., GIFT_db) %>%
    as.data.frame()

community_elements <- predY %>%
  group_by(dominance, genome) %>%
  mutate(row_id = row_number()) %>%
  pivot_wider(names_from = genome, values_from = value) %>%
  ungroup() %>%
  group_split(row_id) %>%
  as.list() %>%
  lapply(., FUN = function(x){x %>%
    select(-row_id) %>%
    column_to_rownames(var = "dominance") %>%
    as.data.frame() %>%
    exp() %>%
    t() %>%
    tss() %>%
    to.community(elements_table,.,GIFT_db) %>% 
    as.data.frame() %>%
    rownames_to_column(var="dominance")
   })

calculate_slope <- function(x) {
  lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
  coef(lm_fit)[2]
}

element_predictions <- map_dfc(community_elements, function(mat) {
      mat %>%
        column_to_rownames(var = "dominance") %>%
        t() %>%
        as.data.frame() %>%
        rowwise() %>%
        mutate(slope = calculate_slope(c_across(everything()))) %>%
        select(slope) }) %>%
      t() %>%
      as.data.frame() %>%
      set_names(colnames(community_elements[[1]])[-1]) %>%
      rownames_to_column(var="iteration") %>%
      pivot_longer(!iteration, names_to="trait",values_to="value") %>%
      group_by(trait) %>%
      summarise(mean=mean(value),
        p1 = quantile(value, probs = 0.1),
        p9 = quantile(value, probs = 0.9)) %>%
      arrange(-mean)

element_predictions %>%
      tt()

trait	mean	p1	p9
B0706	3.719394e-03	-1.082073e-03	9.783638e-03
D0511	3.563852e-03	-9.933616e-04	9.487401e-03
B0708	2.704763e-03	-1.214157e-03	7.062100e-03
D0512	2.402681e-03	-1.952187e-03	7.280814e-03
D0602	2.373480e-03	-1.783097e-03	6.974937e-03
B0214	2.347823e-03	-2.081313e-03	7.224168e-03
D0910	2.231951e-03	-2.563086e-03	7.593477e-03
D0613	2.147485e-03	-4.563149e-03	8.418655e-03
D0203	2.004302e-03	-9.760456e-04	7.434638e-03
D0305	1.910159e-03	-4.180024e-04	6.100309e-03
D0207	1.809096e-03	-2.725130e-03	9.868014e-03
D0513	1.801329e-03	-1.544067e-03	5.819561e-03
B0604	1.617647e-03	-4.135655e-03	1.037624e-02
B0802	1.378116e-03	-1.246028e-03	4.246333e-03
D0308	1.363584e-03	-3.558773e-03	6.831814e-03
B0705	1.111302e-03	-8.095272e-04	3.515149e-03
B0102	1.095401e-03	-8.663080e-04	3.674662e-03
B0704	1.076907e-03	-3.582996e-03	5.682505e-03
B0401	1.044166e-03	-1.642539e-03	5.367512e-03
B0805	1.043874e-03	-6.350291e-04	3.900947e-03
D0309	1.019947e-03	-1.709968e-03	4.691992e-03
D0210	9.688237e-04	-2.354167e-03	4.933620e-03
D0209	9.180974e-04	-2.284493e-03	4.192509e-03
B0702	9.132865e-04	-2.642259e-03	4.471595e-03
B0104	8.786755e-04	-2.579840e-03	4.124706e-03
B0106	8.752008e-04	-1.420667e-03	3.025808e-03
D0706	8.338948e-04	-1.315212e-03	3.192421e-03
D0206	7.887896e-04	-4.538717e-03	8.392202e-03
D0205	7.765677e-04	-1.934373e-03	4.066164e-03
B0707	7.734645e-04	-9.513393e-04	2.493948e-03
D0702	7.323596e-04	-2.619426e-03	4.553279e-03
D0501	6.849191e-04	-2.020073e-03	3.209962e-03
B0803	6.548590e-04	-4.012476e-04	1.464784e-03
B0101	5.800704e-04	-1.155094e-03	2.471310e-03
B0103	5.592310e-04	-3.361862e-05	8.287607e-04
D0213	5.511067e-04	-2.515531e-03	5.318496e-03
D0302	5.458030e-04	-1.781191e-03	3.115639e-03
B0701	5.295642e-04	-3.916704e-03	4.747361e-03
D0306	4.992974e-04	-3.064801e-03	4.561530e-03
D0705	4.872583e-04	-1.271811e-03	3.094388e-03
D0212	4.776537e-04	-2.736532e-03	4.077202e-03
B0703	4.368363e-04	-1.029810e-03	1.933889e-03
B0207	4.305318e-04	-2.461688e-03	3.439118e-03
B0402	3.941157e-04	-2.057514e-03	3.662511e-03
B0403	3.852287e-04	-1.156370e-03	2.864535e-03
S0105	3.690705e-04	-5.329064e-03	4.846884e-03
B0105	3.295090e-04	-1.922728e-03	2.171346e-03
B0216	3.164157e-04	-2.377783e-03	2.971891e-03
D0509	2.828560e-04	-1.800010e-03	2.011564e-03
D0310	2.627430e-04	-4.508906e-03	5.566014e-03
D0505	2.586766e-04	-4.262960e-04	9.193711e-04
B0711	1.513645e-04	-3.176979e-03	3.021999e-03
D0911	1.302221e-04	-3.865237e-03	4.438318e-03
D0708	1.274526e-04	-1.040834e-03	1.261305e-03
D0905	1.273723e-04	-3.903579e-03	4.371356e-03
D0208	1.011182e-04	-4.842200e-04	7.156132e-04
B0712	9.831024e-05	-6.081542e-04	7.699420e-04
B0601	5.982056e-05	-4.426407e-03	3.501262e-03
D0609	3.237726e-05	-1.696310e-03	1.974408e-03
D0304	2.669350e-05	-2.128722e-03	2.169057e-03
D0516	2.131797e-05	-1.434246e-04	2.009280e-04
S0103	1.301808e-05	-1.340765e-03	1.892176e-03
D0801	6.450918e-06	-4.152279e-05	5.871871e-05
D0802	6.450918e-06	-4.152279e-05	5.871871e-05
D0506	4.756104e-06	-1.259604e-03	1.171744e-03
B0310	0.000000e+00	0.000000e+00	0.000000e+00
B1006	0.000000e+00	0.000000e+00	0.000000e+00
B1008	0.000000e+00	0.000000e+00	0.000000e+00
B1011	0.000000e+00	0.000000e+00	0.000000e+00
B1014	0.000000e+00	0.000000e+00	0.000000e+00
B1021	0.000000e+00	0.000000e+00	0.000000e+00
B1022	0.000000e+00	0.000000e+00	0.000000e+00
B1024	0.000000e+00	0.000000e+00	0.000000e+00
B1026	0.000000e+00	0.000000e+00	0.000000e+00
B1029	0.000000e+00	0.000000e+00	0.000000e+00
B1041	0.000000e+00	0.000000e+00	0.000000e+00
B1042	0.000000e+00	0.000000e+00	0.000000e+00
D0303	0.000000e+00	0.000000e+00	0.000000e+00
D0401	0.000000e+00	0.000000e+00	0.000000e+00
D0402	0.000000e+00	0.000000e+00	0.000000e+00
D0403	0.000000e+00	0.000000e+00	0.000000e+00
D0404	0.000000e+00	0.000000e+00	0.000000e+00
D0405	0.000000e+00	0.000000e+00	0.000000e+00
D0406	0.000000e+00	0.000000e+00	0.000000e+00
D0407	0.000000e+00	0.000000e+00	0.000000e+00
D0408	0.000000e+00	0.000000e+00	0.000000e+00
D0612	0.000000e+00	0.000000e+00	0.000000e+00
D0709	0.000000e+00	0.000000e+00	0.000000e+00
D0804	0.000000e+00	0.000000e+00	0.000000e+00
D0808	0.000000e+00	0.000000e+00	0.000000e+00
D0809	0.000000e+00	0.000000e+00	0.000000e+00
D0810	0.000000e+00	0.000000e+00	0.000000e+00
D0812	0.000000e+00	0.000000e+00	0.000000e+00
D0813	0.000000e+00	0.000000e+00	0.000000e+00
D0814	0.000000e+00	0.000000e+00	0.000000e+00
D0815	0.000000e+00	0.000000e+00	0.000000e+00
D0818	0.000000e+00	0.000000e+00	0.000000e+00
D0819	-1.128524e-08	-3.256512e-08	5.342258e-09
D0806	-1.786977e-07	-4.247324e-07	2.379716e-09
B0902	-3.936712e-07	-9.719118e-07	1.588513e-10
D0805	-5.415083e-07	-1.287068e-06	7.211261e-09
D0904	-5.415083e-07	-1.287068e-06	7.211261e-09
D0906	-1.001079e-06	-5.171304e-06	2.727295e-06
S0101	-1.186867e-05	-1.277514e-03	1.079687e-03
D0608	-1.824636e-05	-1.205095e-04	8.309584e-05
D0603	-1.872992e-05	-5.102013e-05	3.486669e-07
B0309	-2.412538e-05	-5.714720e-05	4.375137e-06
D0517	-2.785989e-05	-7.717379e-05	1.480476e-06
B0901	-2.899953e-05	-1.028402e-04	4.716877e-05
B1004	-3.019169e-05	-4.592795e-04	3.819057e-04
B0801	-3.173693e-05	-8.502887e-05	1.035467e-06
D0701	-3.876392e-05	-1.428588e-04	5.790013e-05
D0301	-4.700409e-05	-1.562630e-04	7.031765e-05
D0510	-5.488435e-05	-1.966205e-04	6.508789e-05
D0611	-5.511640e-05	-1.508825e-04	1.576515e-06
B0219	-5.519344e-05	-1.511421e-04	1.447406e-06
D0903	-5.565791e-05	-1.524088e-04	1.443728e-06
B0302	-6.085094e-05	-1.633548e-04	1.617356e-05
B0709	-6.456286e-05	-2.973845e-04	1.815925e-04
D0504	-7.957301e-05	-4.498406e-04	2.593235e-04
B0603	-9.797564e-05	-2.223857e-03	1.778317e-03
S0202	-1.019223e-04	-5.213934e-04	3.112737e-04
D0507	-1.051093e-04	-7.210777e-04	4.531470e-04
B0903	-1.243914e-04	-2.973567e-04	3.768582e-05
D0508	-1.462072e-04	-3.503256e-04	4.595512e-05
D0211	-1.710191e-04	-1.694610e-03	1.521611e-03
D0606	-1.733665e-04	-5.127945e-04	1.671022e-04
D0811	-1.840342e-04	-7.083402e-04	4.729686e-07
D0518	-1.845947e-04	-5.571461e-04	2.020572e-04
D0610	-2.205385e-04	-9.523502e-04	5.196416e-04
B1012	-2.255247e-04	-8.822497e-04	4.519514e-04
D0102	-2.432509e-04	-1.834852e-03	1.151828e-03
D0807	-2.438626e-04	-7.590952e-04	1.652275e-04
D0817	-2.766787e-04	-7.886744e-04	2.460083e-04
B0220	-3.179502e-04	-1.120941e-03	4.809220e-04
D0204	-3.187804e-04	-3.598800e-03	2.226626e-03
D0103	-3.943390e-04	-9.333814e-04	3.989209e-05
B1028	-4.516063e-04	-1.636304e-03	1.024123e-04
D0912	-5.046648e-04	-1.676752e-03	2.640423e-04
S0104	-5.117795e-04	-3.202793e-03	2.649003e-03
B0804	-5.460138e-04	-3.078802e-03	1.486326e-03
B0602	-5.638494e-04	-3.010411e-03	1.127969e-03
D0816	-5.768272e-04	-1.684168e-03	6.426880e-04
B0303	-6.072869e-04	-2.127889e-03	2.104766e-04
B0710	-6.085885e-04	-3.735007e-03	1.240995e-03
D0202	-6.363941e-04	-2.997024e-03	9.871894e-04
D0307	-6.460107e-04	-2.642114e-03	1.292023e-03
B0605	-6.848580e-04	-5.367350e-03	3.765253e-03
B0204	-7.057948e-04	-3.131555e-03	1.329580e-03
D0901	-7.332066e-04	-8.357783e-03	5.384527e-03
B0211	-7.741293e-04	-3.446581e-03	1.173415e-03
D0604	-8.167732e-04	-6.884440e-03	5.889754e-03
D0607	-8.995075e-04	-2.924639e-03	-4.708592e-05
S0301	-9.067269e-04	-5.576003e-03	3.987113e-03
S0201	-1.009194e-03	-5.548315e-03	4.329901e-03
D0104	-1.041557e-03	-2.457207e-03	-2.798915e-06
D0503	-1.567236e-03	-5.906522e-03	-8.936803e-06
D0908	-1.579890e-03	-5.577276e-03	2.538351e-03
B0205	-1.678581e-03	-5.783413e-03	1.406489e-03
B0206	-1.754403e-03	-5.453745e-03	1.397877e-03
D0201	-1.798554e-03	-7.373820e-03	9.412366e-04
D0101	-1.809140e-03	-5.273329e-03	3.698673e-04
D0502	-1.939703e-03	-7.112041e-03	2.495029e-03
B0208	-2.113857e-03	-6.722737e-03	2.027941e-03
B0209	-2.183883e-03	-6.551978e-03	1.810695e-03
B0307	-2.312570e-03	-5.486055e-03	1.613900e-04
B0215	-2.397289e-03	-7.768000e-03	2.276592e-03
B0217	-2.436759e-03	-7.676018e-03	5.391933e-04
B0213	-2.488866e-03	-7.692649e-03	1.778536e-03
D0704	-2.615551e-03	-8.486583e-03	2.006378e-03
B0218	-2.927353e-03	-8.807174e-03	-1.088530e-05
B0210	-2.964770e-03	-8.847444e-03	2.145113e-03
B0212	-3.039284e-03	-7.964048e-03	1.004834e-03
D0902	-3.135033e-03	-1.172783e-02	1.181506e-05
D0601	-3.138992e-03	-1.177955e-02	-1.176429e-05
B0221	-3.948251e-03	-9.488139e-03	2.120444e-04
D0907	-4.549958e-03	-1.067803e-02	-4.870143e-07

# Positively associated
element_predictions %>%
  arrange(-mean) %>%
  filter(mean >0) %>%
  tt()

trait	mean	p1	p9
B0706	3.719394e-03	-1.082073e-03	9.783638e-03
D0511	3.563852e-03	-9.933616e-04	9.487401e-03
B0708	2.704763e-03	-1.214157e-03	7.062100e-03
D0512	2.402681e-03	-1.952187e-03	7.280814e-03
D0602	2.373480e-03	-1.783097e-03	6.974937e-03
B0214	2.347823e-03	-2.081313e-03	7.224168e-03
D0910	2.231951e-03	-2.563086e-03	7.593477e-03
D0613	2.147485e-03	-4.563149e-03	8.418655e-03
D0203	2.004302e-03	-9.760456e-04	7.434638e-03
D0305	1.910159e-03	-4.180024e-04	6.100309e-03
D0207	1.809096e-03	-2.725130e-03	9.868014e-03
D0513	1.801329e-03	-1.544067e-03	5.819561e-03
B0604	1.617647e-03	-4.135655e-03	1.037624e-02
B0802	1.378116e-03	-1.246028e-03	4.246333e-03
D0308	1.363584e-03	-3.558773e-03	6.831814e-03
B0705	1.111302e-03	-8.095272e-04	3.515149e-03
B0102	1.095401e-03	-8.663080e-04	3.674662e-03
B0704	1.076907e-03	-3.582996e-03	5.682505e-03
B0401	1.044166e-03	-1.642539e-03	5.367512e-03
B0805	1.043874e-03	-6.350291e-04	3.900947e-03
D0309	1.019947e-03	-1.709968e-03	4.691992e-03
D0210	9.688237e-04	-2.354167e-03	4.933620e-03
D0209	9.180974e-04	-2.284493e-03	4.192509e-03
B0702	9.132865e-04	-2.642259e-03	4.471595e-03
B0104	8.786755e-04	-2.579840e-03	4.124706e-03
B0106	8.752008e-04	-1.420667e-03	3.025808e-03
D0706	8.338948e-04	-1.315212e-03	3.192421e-03
D0206	7.887896e-04	-4.538717e-03	8.392202e-03
D0205	7.765677e-04	-1.934373e-03	4.066164e-03
B0707	7.734645e-04	-9.513393e-04	2.493948e-03
D0702	7.323596e-04	-2.619426e-03	4.553279e-03
D0501	6.849191e-04	-2.020073e-03	3.209962e-03
B0803	6.548590e-04	-4.012476e-04	1.464784e-03
B0101	5.800704e-04	-1.155094e-03	2.471310e-03
B0103	5.592310e-04	-3.361862e-05	8.287607e-04
D0213	5.511067e-04	-2.515531e-03	5.318496e-03
D0302	5.458030e-04	-1.781191e-03	3.115639e-03
B0701	5.295642e-04	-3.916704e-03	4.747361e-03
D0306	4.992974e-04	-3.064801e-03	4.561530e-03
D0705	4.872583e-04	-1.271811e-03	3.094388e-03
D0212	4.776537e-04	-2.736532e-03	4.077202e-03
B0703	4.368363e-04	-1.029810e-03	1.933889e-03
B0207	4.305318e-04	-2.461688e-03	3.439118e-03
B0402	3.941157e-04	-2.057514e-03	3.662511e-03
B0403	3.852287e-04	-1.156370e-03	2.864535e-03
S0105	3.690705e-04	-5.329064e-03	4.846884e-03
B0105	3.295090e-04	-1.922728e-03	2.171346e-03
B0216	3.164157e-04	-2.377783e-03	2.971891e-03
D0509	2.828560e-04	-1.800010e-03	2.011564e-03
D0310	2.627430e-04	-4.508906e-03	5.566014e-03
D0505	2.586766e-04	-4.262960e-04	9.193711e-04
B0711	1.513645e-04	-3.176979e-03	3.021999e-03
D0911	1.302221e-04	-3.865237e-03	4.438318e-03
D0708	1.274526e-04	-1.040834e-03	1.261305e-03
D0905	1.273723e-04	-3.903579e-03	4.371356e-03
D0208	1.011182e-04	-4.842200e-04	7.156132e-04
B0712	9.831024e-05	-6.081542e-04	7.699420e-04
B0601	5.982056e-05	-4.426407e-03	3.501262e-03
D0609	3.237726e-05	-1.696310e-03	1.974408e-03
D0304	2.669350e-05	-2.128722e-03	2.169057e-03
D0516	2.131797e-05	-1.434246e-04	2.009280e-04
S0103	1.301808e-05	-1.340765e-03	1.892176e-03
D0801	6.450918e-06	-4.152279e-05	5.871871e-05
D0802	6.450918e-06	-4.152279e-05	5.871871e-05
D0506	4.756104e-06	-1.259604e-03	1.171744e-03

# Negatively associated
element_predictions %>%
  arrange(mean) %>%
  filter(mean <0) %>%
  tt()

trait	mean	p1	p9
D0907	-4.549958e-03	-1.067803e-02	-4.870143e-07
B0221	-3.948251e-03	-9.488139e-03	2.120444e-04
D0601	-3.138992e-03	-1.177955e-02	-1.176429e-05
D0902	-3.135033e-03	-1.172783e-02	1.181506e-05
B0212	-3.039284e-03	-7.964048e-03	1.004834e-03
B0210	-2.964770e-03	-8.847444e-03	2.145113e-03
B0218	-2.927353e-03	-8.807174e-03	-1.088530e-05
D0704	-2.615551e-03	-8.486583e-03	2.006378e-03
B0213	-2.488866e-03	-7.692649e-03	1.778536e-03
B0217	-2.436759e-03	-7.676018e-03	5.391933e-04
B0215	-2.397289e-03	-7.768000e-03	2.276592e-03
B0307	-2.312570e-03	-5.486055e-03	1.613900e-04
B0209	-2.183883e-03	-6.551978e-03	1.810695e-03
B0208	-2.113857e-03	-6.722737e-03	2.027941e-03
D0502	-1.939703e-03	-7.112041e-03	2.495029e-03
D0101	-1.809140e-03	-5.273329e-03	3.698673e-04
D0201	-1.798554e-03	-7.373820e-03	9.412366e-04
B0206	-1.754403e-03	-5.453745e-03	1.397877e-03
B0205	-1.678581e-03	-5.783413e-03	1.406489e-03
D0908	-1.579890e-03	-5.577276e-03	2.538351e-03
D0503	-1.567236e-03	-5.906522e-03	-8.936803e-06
D0104	-1.041557e-03	-2.457207e-03	-2.798915e-06
S0201	-1.009194e-03	-5.548315e-03	4.329901e-03
S0301	-9.067269e-04	-5.576003e-03	3.987113e-03
D0607	-8.995075e-04	-2.924639e-03	-4.708592e-05
D0604	-8.167732e-04	-6.884440e-03	5.889754e-03
B0211	-7.741293e-04	-3.446581e-03	1.173415e-03
D0901	-7.332066e-04	-8.357783e-03	5.384527e-03
B0204	-7.057948e-04	-3.131555e-03	1.329580e-03
B0605	-6.848580e-04	-5.367350e-03	3.765253e-03
D0307	-6.460107e-04	-2.642114e-03	1.292023e-03
D0202	-6.363941e-04	-2.997024e-03	9.871894e-04
B0710	-6.085885e-04	-3.735007e-03	1.240995e-03
B0303	-6.072869e-04	-2.127889e-03	2.104766e-04
D0816	-5.768272e-04	-1.684168e-03	6.426880e-04
B0602	-5.638494e-04	-3.010411e-03	1.127969e-03
B0804	-5.460138e-04	-3.078802e-03	1.486326e-03
S0104	-5.117795e-04	-3.202793e-03	2.649003e-03
D0912	-5.046648e-04	-1.676752e-03	2.640423e-04
B1028	-4.516063e-04	-1.636304e-03	1.024123e-04
D0103	-3.943390e-04	-9.333814e-04	3.989209e-05
D0204	-3.187804e-04	-3.598800e-03	2.226626e-03
B0220	-3.179502e-04	-1.120941e-03	4.809220e-04
D0817	-2.766787e-04	-7.886744e-04	2.460083e-04
D0807	-2.438626e-04	-7.590952e-04	1.652275e-04
D0102	-2.432509e-04	-1.834852e-03	1.151828e-03
B1012	-2.255247e-04	-8.822497e-04	4.519514e-04
D0610	-2.205385e-04	-9.523502e-04	5.196416e-04
D0518	-1.845947e-04	-5.571461e-04	2.020572e-04
D0811	-1.840342e-04	-7.083402e-04	4.729686e-07
D0606	-1.733665e-04	-5.127945e-04	1.671022e-04
D0211	-1.710191e-04	-1.694610e-03	1.521611e-03
D0508	-1.462072e-04	-3.503256e-04	4.595512e-05
B0903	-1.243914e-04	-2.973567e-04	3.768582e-05
D0507	-1.051093e-04	-7.210777e-04	4.531470e-04
S0202	-1.019223e-04	-5.213934e-04	3.112737e-04
B0603	-9.797564e-05	-2.223857e-03	1.778317e-03
D0504	-7.957301e-05	-4.498406e-04	2.593235e-04
B0709	-6.456286e-05	-2.973845e-04	1.815925e-04
B0302	-6.085094e-05	-1.633548e-04	1.617356e-05
D0903	-5.565791e-05	-1.524088e-04	1.443728e-06
B0219	-5.519344e-05	-1.511421e-04	1.447406e-06
D0611	-5.511640e-05	-1.508825e-04	1.576515e-06
D0510	-5.488435e-05	-1.966205e-04	6.508789e-05
D0301	-4.700409e-05	-1.562630e-04	7.031765e-05
D0701	-3.876392e-05	-1.428588e-04	5.790013e-05
B0801	-3.173693e-05	-8.502887e-05	1.035467e-06
B1004	-3.019169e-05	-4.592795e-04	3.819057e-04
B0901	-2.899953e-05	-1.028402e-04	4.716877e-05
D0517	-2.785989e-05	-7.717379e-05	1.480476e-06
B0309	-2.412538e-05	-5.714720e-05	4.375137e-06
D0603	-1.872992e-05	-5.102013e-05	3.486669e-07
D0608	-1.824636e-05	-1.205095e-04	8.309584e-05
S0101	-1.186867e-05	-1.277514e-03	1.079687e-03
D0906	-1.001079e-06	-5.171304e-06	2.727295e-06
D0805	-5.415083e-07	-1.287068e-06	7.211261e-09
D0904	-5.415083e-07	-1.287068e-06	7.211261e-09
B0902	-3.936712e-07	-9.719118e-07	1.588513e-10
D0806	-1.786977e-07	-4.247324e-07	2.379716e-09
D0819	-1.128524e-08	-3.256512e-08	5.342258e-09

community_elements %>%
    bind_rows() %>%
    pivot_longer(-dominance, names_to = "trait", values_to = "value") %>%
  filter(trait %in% element_predictions$trait) %>%
  mutate(trait=factor(trait, levels=element_predictions$trait)) %>%
    mutate(dominance=as.numeric(dominance)) %>%
    ggplot(aes(x=dominance, y=value)) +
          geom_smooth(method = lm, formula = y ~ splines::bs(x, 3), se = TRUE) +
          facet_wrap(~trait, ncol=6, scales="free") +
          theme_minimal() +
          labs(x="Dominance",y="Metabolic Capacity Index")

11.2.2.2 Function level

functions_table <- elements_table %>%
    to.functions(., GIFT_db) %>%
    as.data.frame()

community_functions <- predY %>%
  group_by(dominance, genome) %>%
  mutate(row_id = row_number()) %>%
  pivot_wider(names_from = genome, values_from = value) %>%
  ungroup() %>%
  group_split(row_id) %>%
  as.list() %>%
  lapply(., FUN = function(x){x %>%
    select(-row_id) %>%
    column_to_rownames(var = "dominance") %>%
    as.data.frame() %>%
    exp() %>%
    t() %>%
    tss() %>%
    to.community(functions_table,.,GIFT_db) %>% 
    as.data.frame() %>%
    rownames_to_column(var="dominance")
   })

#max-min option
calculate_slope <- function(x) {
  lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
  coef(lm_fit)[2]
}

function_predictions <- map_dfc(community_functions, function(mat) {
      mat %>%
        column_to_rownames(var = "dominance") %>%
        t() %>%
        as.data.frame() %>%
        rowwise() %>%
        mutate(slope = calculate_slope(c_across(everything()))) %>%
        select(slope) }) %>%
      t() %>%
      as.data.frame() %>%
      set_names(colnames(community_functions[[1]])[-1]) %>%
      rownames_to_column(var="iteration") %>%
      pivot_longer(!iteration, names_to="trait",values_to="value") %>%
      group_by(trait) %>%
      summarise(mean=mean(value),
        p1 = quantile(value, probs = 0.1),
        p9 = quantile(value, probs = 0.9)) %>%
      arrange(-mean)

function_predictions %>%
      tt()

trait	mean	p1	p9
B07	1.244660e-03	-0.0006954170	3.339535e-03
B08	7.258465e-04	-0.0001971510	1.856081e-03
B01	7.196813e-04	-0.0009349156	2.121851e-03
B04	5.722166e-04	-0.0015697967	3.880385e-03
D05	5.607811e-04	-0.0004522992	1.660790e-03
D03	5.606577e-04	-0.0012031156	2.944461e-03
D02	4.208314e-04	-0.0014837961	2.708039e-03
S01	1.046572e-04	-0.0022765478	2.167428e-03
D04	0.000000e+00	0.0000000000	0.000000e+00
B09	-3.445384e-05	-0.0000976261	3.119699e-05
B10	-5.518709e-05	-0.0001494436	8.242006e-06
D08	-7.177074e-05	-0.0001820432	2.008773e-05
B06	-1.850995e-04	-0.0021286727	1.386341e-03
D06	-3.221232e-04	-0.0018931867	7.764011e-04
S02	-4.648308e-04	-0.0024154945	1.803723e-03
B03	-5.109473e-04	-0.0011883082	1.953477e-06
D07	-6.897856e-04	-0.0026307024	1.062556e-03
D09	-7.336735e-04	-0.0022399454	6.072017e-04
S03	-9.067269e-04	-0.0055760032	3.987113e-03
D01	-1.059485e-03	-0.0028771127	1.623250e-05
B02	-1.502756e-03	-0.0043474039	7.369526e-04

# Positively associated
function_predictions %>%
  arrange(-mean) %>%
  filter(mean >0) %>%
  tt()

trait	mean	p1	p9
B07	0.0012446598	-0.0006954170	0.003339535
B08	0.0007258465	-0.0001971510	0.001856081
B01	0.0007196813	-0.0009349156	0.002121851
B04	0.0005722166	-0.0015697967	0.003880385
D05	0.0005607811	-0.0004522992	0.001660790
D03	0.0005606577	-0.0012031156	0.002944461
D02	0.0004208314	-0.0014837961	0.002708039
S01	0.0001046572	-0.0022765478	0.002167428

# Negatively associated
function_predictions %>%
  arrange(mean) %>%
  filter(mean <0) %>%
  tt()

trait	mean	p1	p9
B02	-1.502756e-03	-0.0043474039	7.369526e-04
D01	-1.059485e-03	-0.0028771127	1.623250e-05
S03	-9.067269e-04	-0.0055760032	3.987113e-03
D09	-7.336735e-04	-0.0022399454	6.072017e-04
D07	-6.897856e-04	-0.0026307024	1.062556e-03
B03	-5.109473e-04	-0.0011883082	1.953477e-06
S02	-4.648308e-04	-0.0024154945	1.803723e-03
D06	-3.221232e-04	-0.0018931867	7.764011e-04
B06	-1.850995e-04	-0.0021286727	1.386341e-03
D08	-7.177074e-05	-0.0001820432	2.008773e-05
B10	-5.518709e-05	-0.0001494436	8.242006e-06
B09	-3.445384e-05	-0.0000976261	3.119699e-05

#lm option (incoherent with Bayesian framework)
#function_predictions <- community_functions %>%
#   bind_rows() %>%
#   pivot_longer(-dominance, names_to = "trait", values_to = "value") %>%
#  group_by(trait) %>%
#  mutate(dominance=as.numeric(dominance)) %>%
#  summarize(model=list(lm(value ~ dominance))) %>%
#  ungroup() %>%
#  select(trait,model) %>%
#  mutate(estimate = map_dbl(model, ~broom::tidy(.) %>% filter(term == "dominance") %>% pull(estimate))) %>%
#  mutate(p_value = map_dbl(model, ~broom::tidy(.) %>% filter(term == "dominance") %>% pull(p.value))) %>%
#  mutate(p_value_adj = p.adjust(p_value, method = "bonferroni")) %>%
#  select(-model) %>%
#  arrange(-estimate) %>%
#  filter(p_value_adj < 0.05) 

#brms option (very lengthy) > even more if group is added as random effect (not in this example)
#function_predictions <- community_functions %>%
#   bind_rows() %>%
#   pivot_longer(-dominance, names_to = "trait", values_to = "value") %>%
#  mutate(dominance=as.numeric(dominance)) %>%
#  group_by(trait) %>%
#  do(model=list(brm(value ~ dominance, data=.))) %>%
#  ungroup() %>%
#  select(trait,model) %>%
#  mutate(estimate = map_dbl(model, ~broom.mixed::tidy(.) %>% filter(term == "dominance") %>% pull(estimate))) %>%
#  mutate(error = map_dbl(model, ~broom.mixed::tidy(.) %>% filter(term == "dominance") %>% pull(std.error))) %>%
#  mutate(error = map_dbl(model, ~broom.mixed::tidy(.) %>% filter(term == "dominance") %>% pull(conf.low))) %>%
#  mutate(error = map_dbl(model, ~broom.mixed::tidy(.) %>% filter(term == "dominance") %>% pull(conf.high))) %>%
#  arrange(-estimate)

community_functions %>%
    bind_rows() %>%
    pivot_longer(-dominance, names_to = "trait", values_to = "value") %>%
  filter(trait %in% function_predictions$trait) %>%
  mutate(trait=factor(trait, levels=function_predictions$trait)) %>%
    mutate(dominance=as.numeric(dominance)) %>%
    ggplot(aes(x=dominance, y=value)) +
          geom_smooth(method = lm, formula = y ~ x, se = TRUE) +
          #geom_smooth(method = lm, formula = y ~ splines::bs(x, 3), se = TRUE) +
          facet_wrap(~trait, ncol=5, scales="free") +
          theme_minimal() +
          labs(x="Dominance",y="Metabolic Capacity Index")